
function boole = system_has_solution(dataset,epsilons,totobs,nobs,ignore_errors)

	if nargin<5
		ignore_errors=0;
	end
	NN = -eye(4*nobs); b_nn = -0.0001*ones(4*nobs,1); % non negativity constraints

	for i = 1:size(dataset,1)

		epsilon = epsilons(i);

		xc = table2array(dataset(i,1:totobs))';
		yc = table2array(dataset(i,totobs+1:2*totobs))';
		px = table2array(dataset(i,2*totobs+1:3*totobs))';
		py = table2array(dataset(i,3*totobs+1:4*totobs))';

		D = [xc yc px py];

		% equality constraints linked to the fact that i choose the same value in M at different decisions
		M = [xc(1:nobs); yc(1:nobs)];
	    m = M; k = 0;
	    Veq = zeros(1,2*nobs);
	    while ~isempty(m)>0
	    	same_idx = find(M==m(1));
	    	m(m==M(same_idx(1))) = [];
	    	
	    	for j=2:length(same_idx)  % note: this runs only if length(same_idx)>=2
	    		k = k+1;
	    		Veq(k,[same_idx(1) same_idx(j)]) = [1 -1];
	        end
	    end
	    Veq = [zeros(size(Veq,1),2*nobs) Veq];
	    b_veq = zeros(size(Veq,1),1);

		A1 = [(1-epsilon)*eye(nobs) -eye(nobs) zeros(nobs) zeros(nobs);
			   -eye(nobs) (1-epsilon)*eye(nobs) zeros(nobs) zeros(nobs)];

		b1 = zeros(2*nobs,1);

		[A2, A2eq, b2, b2eq] = write_constraints_part3(D,1,1,nobs);
		[A3, A3eq, b3, b3eq] = write_constraints_part3(D,1,2,nobs);
		[A4, A4eq, b4, b4eq] = write_constraints_part3(D,2,1,nobs);
		[A5, A5eq, b5, b5eq] = write_constraints_part3(D,2,2,nobs);

		A = [NN;A1;A2;A3;A4;A5] ; b = [b_nn;b1;b2;b3;b4;b5];
		Aeq = [Veq;A2eq;A3eq;A4eq;A5eq]; beq = [b_veq;b2eq;b3eq;b4eq;b5eq];

	    options = optimoptions('linprog','Algorithm','interior-point');
	    [~, ~, exitflag] = linprog(zeros(4*nobs,1),A,b,Aeq,beq,[],[],options);

		if exitflag == 1
			boole(:,i) = 1;
		elseif exitflag == -2
			boole(:,i) = 0;
        else
        	if ignore_errors
            	boole(:,i) = -1;
            else
            	error("some other problem (e.g. didn't find a solution but one might exist) \n code: %d \n individual %d", exitflag, i);
			end
		end
	end
end

function [A, Aeq, b, beq] =  write_constraints_part3(D,l,l_prime,nobs)

	A = zeros(nobs^2,4*nobs); b = zeros(nobs^2,1); idx_eq = zeros(nobs^2,1);
	n=0;
	for i = 1:nobs
	    for j = 1:nobs
	        n = n+1;
	        if D(j,l_prime)~=D(i,l)
	            A(n, (1+l_prime)*nobs + j ) = 1;
	            A(n, (1+l)*nobs + i ) = A(n, (1+l)*nobs + i ) - 1;
	            A(n, (l-1)*nobs + i ) = -D(i,l+2)*( D(j,l_prime)-D(i,l) );
	            b(n) = -1; % costant = 1
	        else
	        	idx_eq(n) = 1;
	            
	        end
	    end
	end

	Aeq = A(idx_eq==1,:);
	A   = A(idx_eq==0,:);
	beq = b(idx_eq==1);
	b   = b(idx_eq==0);
end
